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Abstract. Silicene consists of a monolayer of silicon atoms in a buckled honeycomb 

structure. It was recently discovered that the symmetry of such a system allows for 
interesting Rashba spin-orbit effects. A perpendicular electric field is able to couple to 
the sublattice pseudospin, making it possible to electrically tune and close the band 
gap. Therefore, external electric fields may generate a topological phase transition 
from a topological insulator to a normal insulator (or semimetal) and vice versa. 
The contribution of the present article to the study of silicene is twofold: First, we 
perform a group theoretical analysis to systematically construct the Hamiltonian in 
the vicinity of the K points of the Brillouin zone and discover a new, but symmetry 
allowed term. Subsequently, we identify a tight binding model that corresponds to the 
group theoretically derived Hamiltonian near the K points. Second, we start from this 
tight binding model to analyze the topological phase diagram of silicene by an explicit 
calculation of the Z2 topological invariant of the band structure. To this end, we 
calculate the Z2 topological invariant of the honeycomb lattice in a manifestly gauge 
invariant way which allows us to include Sz symmetry breaking terms - like Rashba 
spin orbit interaction - into the topological analysis. Interestingly, we find that the 
interplay of two Rashba terms can generate a non-trivial quantum spin Hall phase in 
silicene. This is in sharp contrast to the more extensively studied honeycomb system 
graphene where Rashba spin orbit interaction is known to compete with the quantum 
spin Hall effect in a potentially detrimental way. 
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1. Introduction 

One of the main subjects of current interest in condensed matter physics is the search 
for materials that host topological insulator (TI) phases [1, 2, 3]. Two dimensional TIs 
exhibit the quantum spin Hall effect (QSHE) with gapless edge states and a finite energy 
gap in the bulk [4, 5, 6]. The first proposal of this state of matter was made by Kane 
and Mele [4] on the basis of graphene in the presence of spin-orbit interaction (SOI). 
However, the relevant SOI in graphene turns out to be rather small [7] such that the 
effect seems to be inaccessible in experiments. This situation is different in HgTe/CdTe 
quantum wells where the QSHE was also predicted theoretically [6] and experimentally 
seen soon after [8]. 

Recently, a single layer of silicon atoms - called silicene - has been synthesized 
exhibiting an analogous honeycomb structure as graphene [9, 10, 11]. Since silicon is 
heavier than carbon, the spin-orbit gap in silicene is much larger than in graphene. 
Therefore, if it was possible at some point to prepare clean silicene, it should be feasible 
to experimentally access the QSHE in this material. Similar to graphene, the unit cell 
of silicene contains two atoms which gives rise to two different sublattices. In contrast 
to graphene, however, the silicene sublattices are found to be arranged in a buckled 
structure pointing out-of-plane [12]. Due to the broken sublattice symmetry, the mobile 
electrons in silicene are therefore able to couple differently to an external electric field 
than the ones in graphene. This difference is the reason for new Rashba spin-orbit 
coupling effects that allow for external tuning and closing of the band gap in silicene [13]. 
Consequently, an electrically induced topological quantum phase transition is possible. 
It is natural to ask whether this phase transition can in principle go both ways, i.e., 
whether the electric field can be used to destroy and generate the QSHE. Refs. [4, 14] 
clearly show that a different potential on the two sublattices of a honeycomb lattice leads 
to a transition from a TI to a trivial insulating state. In Ref. [15], some indications have 
been presented that the interplay of two silicene specific Rashba terms can even induce 
the QSHE starting from a trivial insulating band structure in the absence of Rashba 
SOI. 

The quantum spin Hall (QSH) phase is distinguished from a normal insulating 
phase by a bulk Z2 topological invariant [5, 16]. For a minimal model of the QSHE 
in graphene, this invariant has been analytically calculated in a seminal work by Kane 
and Mele [5]. However, the original formulation of the Z12 invariant in terms of Bloch 
functions does not contain a constructive prescription as to its numerical evaluation. 
Subsequent work on the topological properties of the band structure of silicene was 
restricted to the absence of terms breaking the spin 5'^-conservation or to employing the 
bulk-boundary correspondence in silicene nanoribbons [14, 15]. 

Evidently, a full topological analysis of silicene is missing and, as we show 
below, important to clearly identify phenomenological differences between graphene 
and silicene. In this work, we employ Prodan's method [17] to calculate the topological 
invariant without any further symmetry assumptions in a manifestly gauge invariant 
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way to provide a conclusive analysis of the novel features of silicene regarding QSH 
physics. In particular, we establish that, in contrast to graphene, the QSHE can be 
generated by Rashba SOI in silicene. 

The bulk of this article consists of two parts which we keep fairly self contained to 
allow the reader to follow our analysis a la carte. In Sec. 2, we analyze the symmetries 
of the lattice of silicene. This analysis allows us to mathematically construct the low- 
energy Hamiltonian (close to the K points of the Brillouin zone) by means of the 
invariant expansion method with a particular focus on terms involving a perpendicular 
electric field. Thereby, we discover a new, symmetry allowed term of the low energy 
Hamiltonian. Furthermore, a tight-binding calculation is performed to verify the terms 
previously derived from symmetries and to estimate their magnitude. The reader who is 
more interested in quantum spin Hall physics can directly go to Sec. 3, where we study 
the topological properties of the band structure of silicene by explicitly calculating the Z2 
topological invariant in a manifestly gauge invariant way. In this section, the possibility 
of a topological phase transition induced by an external field is carefully examined which 
enables us to correct previously proposed phase diagrams. Finally, we conclude in Sec. 4. 

2. Symmetry based derivation of the Hamiltonian 
2.1. Identification of the lattice symmetry 

Silicene is a monolayer of silicon atoms arranged in a buckled honeycomb lattice (see 
Fig. 1 for a schematic). In contrast to graphene, the two basis atoms of the unit cell 
(called A and B) are separated perpendicular to the atomic plane at a distance 21 with 
/ = 0.23A [14]. As there is no translation symmetry in the out-of-plane direction, the 
material is quasi-two-dimensional. The buckling is quantified by an angle 6 > 90° as 
shown in Fig. 1. 




Figure 1. Schematic of tlie real lattice of silicene. The sublattices (denoted A and B) 
of the honeycomb structure are spatially separated in z-direction. The buckling-angle 
9 is found to be 101.7° in synthesized silicene [12]. 

In Fig. 2, we provide an illustration of the lattice of silicene - in real and reciprocal 
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space. The latter is again hexagonal rotated by 7r/2. Basis vectors defining the unit cell 
are given by 
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in real space and by 
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in reciprocal space. Here, a is the distance between two neighboring silicon atoms, 
which is equal to the lattice constant divided by a factor \/3. In Fig. 2(b), the buckling 
is visualized by the red dots. Moreover, all potential symmetry axes and planes are 
marked as dashed and dotted fines. 





(a) ib) 

Figure 2. Symmetry operations of the reciprocal (a) and the real lattice (b). The 
operations Cn denote rotations by 27r/n around the z-axis; and C'^ refer to rotations 
around the labeled axis, lying within the atomic plane, a describes reflection planes 
spanned by the labeled axis and the z-axis perpendicular to the plane. The primed 
operations cross atomic sites, while the double-primed ones do not. A and B denote 
different sublattices, K and K' inequivalent corner points of the first Brillouin zone. 
The red dots implicate, that in real space the atomic sites are shifted perpendicular to 
the plane. 



In real space witfiout buckling (i.e. tfie grapfiene case), all symmetries sfiown in 
Fig. 2 will be valid witfi tfie additional symmetry of a reflection at tfie plane itself. 
Tfiis situation corresponds to tfie point group D^h- In tfie presence of tfie buckling (i.e. 
tfie silicene case), tfie atomic plane is no longer a symmetry plane. Furtfiermore, tfie 
symmetries Cg, C2, C2, and a" are broken. Tfie only symmetry classes left tfien are C3, 
C2, and a' wfiicfi we can summarize as point group Cs^ witfi an additional symmetry 

Tfie r point in reciprocal space always exfiibits tfie same symmetries as tfie real 
lattice. Note, fiowever, tfiat tfie same reflection planes and axes are denoted differently 
in real and reciprocal space, differing by a prime in our notation, e.g. O^lr) corresponds 
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to C2{k) (see the caption of Fig. 2 for a careful definition of our notation). Importantly, 
operations in both spaces refer to the same coordinate system. This particularly means, 
that in reciprocal space the classes and a' are forbidden. Thus, at the F point, the 
group of the wave vector is characterized by the symmetry operations C3, and a" 
which corresponds to point group D3 with the additional symmetry class a". At the K 
points of the Brillouin zone, the symmetry of the group of the wave vector is further 
reduced to point group D^. 



2.2. Invariant Expansion 

The full knowledge of the lattice symmetries makes it in principle possible to construct a 
low-energy Hamiltonian by expansion around high-symmetry points. This can be done 
with powerful and well-established approaches, for instance, the invariant expansion 
[18, 19, 20]. In undoped silicene, the Fermi level is located at the i^-points of the 
Brillouin zone. Hence, a low-energy Hamiltonian constructed by a symmetry analysis 
near the i^'-points will capture the essential physical properties of the system. Before we 
perform such a symmetry analysis, let us, for completeness, shortly review the general 
theory of the invariant expansion. 

In this approach, the Hamiltonian is constructed to be invariant under all symmetry 
operations of the point group. When spin is included, it is composed block-wise 
according to irreducible representations (IR) of the corresponding double group. Each 
block ?^"^(/C), that may depend on some general tensor components /C, has to reflect 
the symmetry properties of all the IRs F^, that are contained in the product 

K 

where the integer is the multiplicity. The most general ansatz is then 

U^\1C) = J^^^^^^J"'''^ *• (3) 

Here, k runs over all IR F^ included in the product and a"^''^ are material-specific 
constants. The symmetrized tensor operator components IC^'^'^^ * are parameters like 
the wave vector k, as well as the external magnetic B or electric field S. Furthermore, 
X^'^'^^ are the symmetrized matrices we wish to find to construct the Hamiltonian. 
Evidently, A und n are indices numbering the different possible matrix- and tensor- 
components. These indices run from unity to limits that depend on the order of the 
expansion. One set of the k, A, and fi is called an invariant. 

The irreducible tensor components belonging to an IR are composed analogously 
to its eigenfunctions. With the help of projection operators [21], that contain the 
matrix representations of each IR, we can combine tensor components at any order 
and project them onto the appropriate IR of the point group. To construct the basis 
matrices Xi^'^^ , we can use the Wigner-Eckart theorem [22] . Since any point group is a 
subgroup of the full rotation group TZ, the eigenfunctions of any IR of the point group are 
also eigenfunctions of TZ and can be written in terms of spherical harmonics. Angular 
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momentum quantum numbers are assigned to each IR, by defining the axial vector 
components R^, Ry, and Rz to be the real angular momentum eigenstates with quantum 
number / = 1. Each IR can now be classified with angular momentum quantum numbers 
by comparing its eigenfunctions to a table of real spherical harmonics. A symmetrized 
matrix * (that transforms like the IR of dimension and numbering u = 1, . . . Z^,) 
which is contained in the block T* x Tj is then derived by Clebsch-Gordan coefficients 
(CGC) 



Mm 





J 


i I y 









(4) 



where the order of i, j, and k is crucial. In Eq. (4), I is again the multiplicity, that 
shows how often the IR is included in the product F* x F^. For multiplicities larger 
than one, we find linear independent sets of basis matrices. 

2.3. Application to the ir-band Hamiltonian of silicene near the K points 

We now start to perform an invariant expansion of the group D3, that we found to 
be the group of the wave vector at the i^-points in silicene (see Ref. [23] for a general 
discussion of this point group). We are interested in the two-dimensional vr-bands of 
silicene. The orbital part of the wavefunction corresponds to the two-dimensional IR 
F3 of the group D-^. Since we wish to include spin in our symmetry analysis, we have 
to complement this IR by F4, the appropriate double group IR of D^. Consequently, 
we start from the following product of representations F3 x F4 = F4 -|- F5. The resulting 
(4 X 4) Hamiltonian based on the invariant expansion can then be composed in blocks 
of four (2 X 2) matrices 

■U= I ^"^^ '^^^ ) (5) 

The construction by CGCs (as discussed in the previous section) implies that the vr-band 
Hamiltonian is given in the basis (V'a/3, "ipso:, —ipAOi.'ipBPY', with ipA(B) = Rx T^Ry and 
G {t)i}- Each block of the Hamiltonian is thus decomposed into IRs. General 
tensor components of any desired order can be found with the help of projection 
operators from the point group character table. In Table 1, we list components of 
interest up to third order. Identifying combinations of general tensor components of R 
with angular momentum quantum numbers, we can assign such quantum numbers to 
the IRs itself. The basis matrices are then identified from CGCs according to Eq. (4) 
and listed in Table 2. Note that we have applied a unitary basis transformation to 
guarantee a more symmetric basis {ipA^i i'sCi, '^AO^^ipsPY ■ This transformation is done 
for comparison with the graphene case carefully analyzed in Ref. [24]. We use the 
symmetrized matrices and tensor components listed in Tables 1 and 2 to expand the 
Hamiltonian up to first orders in k and the electric field Ez perpendicular to the plane. 
Note that since "H is Hermitian, the coefficients of the diagonal blocks have to be real, 
while those of the off-diagonal blocks can in principle be imaginary. We write 7 and 
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Table 1. Character table and irreducible tensor components of up to third order 
in the tensor components x. x,y,z denote polar vector components and Rx,Ry,Rz 
axial vector components. Here, axial and polar vector components appear on equal 
footing as the symmetry operations of D3 do not provide mirror planes. Therefore, for 
simplicity, combinations with components of R are not listed again in higher orders 
but they are of course present. 
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Table 2. Symmetrized matrices for blocks arising from the product Fa x r| = r4 + r5 
of the double group of D3. Our choice of basis is {ipAl3,tpBCe,'fpACi,4'BP)'^ ■ denote 
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^7' = 7 to include both cases with real coefficients 7, 7'. Then, we obtain for the 
Hamiltonian 

7^44 = Qiil + a2(JzEz + a^{axkx - (Jyky) + aiE^ia^ky + ayk^), (6) 

^55 = + /?2(Ta + p3(TzE, + pAOyE,, (7) 

"^45 = liijk^j. - la^ky) + ij[{Ik.^ - icr^ky) + (8) 
l2{crxkx + o-yky) + «72(f7x.A';,x- + o-yky) + 
^sEzilky + lOzkx) + i^'^Ez{lky + lOzka:) + 
liEz{(Txky - Uykx) + iiiEz{axky - Uykx). 

For a better physical interpretation, we once more change the basis by an unitary 
transformation and present the total Hamiltonian in the basis (■0^^, ■^^a, ■^s^, ■0bq;)-^. 
After this basis transformation, the Pauli matrices a and s obtain the following physical 
meaning: a acts on the space of sublatticcs A and and s on the spin space. In lowest 
order, the Hamiltonian in our new basis exhibits sixteen terms labeled by coefficients oi 
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"H^ = ail + a2(T^Sz + aad^So-E^ + a4(XoS^E^ + a^^a^cSy - (JyS^)Ez + (9) 
a&ic^xSx + OySy) + a-j{axkx + ayky)so + asaois^K + Syky) + 

^°i\'^x{^Sxkx ^yky) (^yiSyk^ + Sxky)] + 

(^10Ez[<^x{Sxky + Sykx) + CTyi^Sxkx Sy/Cy)] + 

auEzCTols^ky - Syk^) + ai2Ez{(yxky - aykx)so + 
O'ls'^zi^xky — Sykx) + ai4(cra;A;y — (Tykx)Sz + 
ai^EzOzisxkx + SyA;^) + aiQEzia^k^ + ayky)sz. 

Note, that the coefficients in Eqs. 6 to 8 and in Eq. 9 are consistent in the sense that 
they are related by mutual linear combinations. Here, we would already like to point 
out that a new term proportional to the coefficient 04 is found, coupling directly the 
out-of-plane components of spin and electric field. 

2.4- Consequences of time-reversal symmetry 

The terms derived so far in the invariant expansion have not yet been checked for 
consistency with time reversal symmetry (TRS) which holds throughout this work since 
we only consider the influence of electric fields. Importantly, the point group D3 + a" 
provides symmetry operations that map the inequivalent corner points K and K' of the 
Brillouin zone onto each other. For example, we can consider the reflection at a plane 
perpendicular to the atomic plane and including the y-axis (called Ry in Ref. [24]). This 
symmetry operation has the following impact on the Hamiltonian 

V{Ry)%^{lC)V{Ry)-^ = n'^'iR-'JC), 

where V{Ry) is the matrix representation of the operation Ry. Polar (x) and axial 
(R) tensor components will then transform like Ry^ : {x,y,z) = {—x,y,z) and 
Ry^ : {Rx, Ry, Rz) = {Rx, ~Ry, —Rz)- The sublattices will not be changed by Ry. In a 
compact notation, we can write Ry = aoSx- With the help of Ry we find the form of the 
Hamiltonian at one of the two inequivalent K points. All terms, that change sign under 
Ry will thus be modified by the additional parameter = ±1 to mark the difference 
between the valleys K and K'. The true time reversal operator equals T = —aoTxSyC, 
where C is the complex conjugation operator [25] and a Pauli matrix corresponding 
to the valley isospin. Its impact on the Hamiltonian can be written as 

where ^ = ±1 depending on whether the tensor component K, changes sign under time 
reversal or not. Hence, both operators Ry and T lead to transformations between the 
valleys. We now combine both of them to a (new) time-reversal operator within a single 
valley 6 (in the spirit of Ref. [24]): Q{Ry) = TVi^Ry) = iSzC. This operator yields an 
additional symmetry constraint for all terms 

en{}C)e-' = SzWm-^ic)sz, (lo) 
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which forces the coefficients ae, as, ag, and ai2 to vanish. Our final result of this section 
is the low-energy Hamiltonian of silicene near the i^'-points (in first orders in k and E^) 
presented in the basis {ipA(3,i'ACi,4'Bf!^,4'BCi)'^ ■ It is given by 

= Oil + a2T^crzSz + a^a^soE^ + a^T^aos^E^ + (11) 
abiTzC^xSy - (TySx)E^ + a^i^T^cr^k^ + ayky)so + 

(^xi^Sxky -\- Sykx) ~\~ Tz(^ yi^Sxkx "^3/^3/)] ~l~ 

aiiE^ao{sxky - Sykx) + ai^a^^Sxky - Sykx) + 

ky)Sz- 

In the above Hamiltonian, the terms proportional to 02, 05, and a-j are well known from 
the graphene literature to describe its fundamental properties near the i^-points [4]. 
Further corrections proportional to aio and an can as well be found in graphene [24]. 
Beyond this, silicene exhibits specific terms proportional to 03 and ai3 as reported in 
Refs. [13, 12]. Additionally, we find a new term proportional to 04 as well as a higher 
order correction to the linear dispersion (proportional to aie). We conclude, that the 
low-energy Hamiltonian of silicene near the i^-points contains interesting SOI terms that 
are absent in graphene because of the difference in the symmetry of the two honeycomb 
lattices. 



2.5. Corresponding tight-binding model 

To supplement the symmetry analysis, a corresponding tight-binding model is discussed 
next. This model enables us to construct a valid Hamiltonian for the full Brillouin zone 
which is crucial for the subsequent topological analysis. 

Let us start with a brief discussion of the internal spin-orbit coupling and 
subsequently introduce the other important terms of the tight-binding model. In real 
space, silicene is described in terms of a lattice model including the (standard) spin-orbit 
coupling term [4, 12] 

where F is the force stemming from the electric potential V , p is the momentum, 
and s the spin of the electron. The (internal) electric force in silicene is provided by 
the crystal field in in-plane and out-of-plane direction. Since the momentum operator 
is oriented along nearest-neighbor or next-nearest neighbor bonds, the silicene lattice 
forbids terms involving a crystal in-plane force coupling to a nearest-neighbor bond. 
This consideration leads to the following spin-orbit Hamiltonian [12] 

2 

Hso = lXso Vijcl^sf Cj (3-1- Xr,2 Yl /^4a(^X C^V? (13) 

with so far undefined parameters Xso and Ar.,2- In the latter equation, a and (3 are 
spin quantum numbers; the indexes i,j label the atomic site/orbital. Here and in 
the following, {i,j) denotes nearest neighbors and next-nearest neighbors. The 
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neighboring sites are each time connected by the vector dij with its corresponding 
unit vector 

4- The sign Vij = ±1 refers to the next-nearest neighbor hopping being 
anticlockwise or clockwise with respect to the positive z-axis. Furthermore, /Xi = ±1 is 
introduced to distinguish between the A{B) site. 

Additionally, there is a regular nearest-neighbor hopping term and on-site energies 
of the form 

Ho = ^^ic\aCia- UA^Cjp, (14) 

where is the on-site energy of the atomic orbital and tij the hopping parameter for 
hopping between the orbitals i and j of neighboring atomic sites. When an external 
electric field is applied perpendicularly to the atomic plane of silicene, a staggered 
sublattice potential of the form 

He = KY^ Eii^icl^da (15) 

ia 

is generated [14] with a so far undefined parameter Ag. Interestingly, El allows for 
on-site transitions between the Pz and s orbitals [26] . 

To describe a full next-nearest neighbor tight-binding model, we also introduce 
spin-orbit coupling terms due to external electric fields, so-called Rashba terms. The 
simplest ones are 

HR = ^Xr,l Yl 4(^"'x4):%i?^+^Ae,2 Yl v.Aasfc,p^i,Ei + 

Y cl{sxd^l^)fc,pEi. (16) 

((ij»;a/3 

Interestingly, the second term proportional to Ag 2 is a multiplicative combination of the 
spin-orbit term and the staggered sublattice potential. This term corresponds to the 
new term that we have found in the invariant expansion model (proportional to a^). 
The full tight-binding Hamiltonian is then given by 

H = Ho + Hso + He + Hr. (17) 

The terms given above will reproduce our Hamiltonian derived by symmetry analysis of 
Eq. (11), when being expanded around the Dirac i^-point. 

2.6. Tight-binding model including it and a-bands of silicene 

To be able to estimate the coupling constants introduced above, we now briefly discuss 
a tight-binding model presented in Ref. [12] including vr and a-bands in silicene using 
the basis 

{\pf), bf), b^), b^), k^), \Py), bf), k^)} X {t,i}- (18) 

Nearest-neighbor hopping matrix elements between given orbitals can then be expressed 
by parameters Vi, V2, and V3 as functions of Slater bond parameters V^p^r, Vppa, Vspa 
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and curvature angle 6 with 
3 
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sin^K 



spai 
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(see Ref. [12] for details of the modeling). Beyond the analysis done in Ref. [12], we 
consider the influence of an electric field applied perpendicularly to the atomic plane. 
A staggered sublattice-Hamiltonian is then induced that takes in the basis (18) the 
following form 
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(19) 



as only transitions between Pz and s orbitals of the same site are allowed [26]. In 
this equation, e is the charge of an electron and zq the Stark element weighting the 
transition. Spin-orbit coupling is now included in the same way as in Ref. [12] by the 
term Hso = ■ s where L is the angular momentum vector, s the spin operator, 
and A the coupling constant. The tight-binding model of vr and a bands allows us to 
estimate some of the parameters introduced in the previous section as a function of the 
parameters Vi_3 (that are in principle known) as well as the Stark element zq. To do 
so, we now apply a unitary transformation U to the Hamiltonian that corresponds to 
the following change of basis 

U^{ptv^.vtvts\.-^ 



z ' Vz , Vy 5 Vx , '-' ' fy , fx , 



X {t,i} = {01,04,02,05,03,06,07,08} X {t,i}, (20) 



with mixed orbitals, for example, 
transformed Hamiltonian 

H' 



UuPz + W2lS^ 



h(Px ~ ))• Then, the 



V2' 



U\Hso + He)U (21) 

can be analyzed perturbatively on the first (4 x 4) block corresponding to the basis 
{01 t, 01 i, 04 t, 04 t}, which is expected to correspond to energy eigenstates near the 
Fermi energy [12]. Listing only terms including the external electric field, we find in 
first and second order perturbation theory the following three terms 

^(1,2) 



XeCTzSo + Ae,2O"0S^ + K,l{(^xSy + CT^Sj 



with parameters 
A. 



2eEzZQUiiU2i 



2eEzZoVs 

Vo 



(22) 



Group theoretical and topological analysis of the quantum spin Hall effect in silicene 12 



X 



e,2 



2Vi 



Ar,l 



AeE^zo 



2V2 

In the last step, the expressions were approximated under the assumption of low 
buckling, meaning 6 — )■ 90° and V3 — > 0, where only the term of lowest order in V3 
was kept. The term proportional to Ag represents the on-site hopping occuring only in a 
buckled structure. With the one proportional to Ae,2, a Rashba term of first order in SOI 
A and is generated that is related to a next-nearest neighbor hopping. It depends 
linearly on V3, so this term is specific for the buckled silicene structure as well. Finally, 
the term proportional to A,.,i is the well-known Rashba spin orbit coupling reported 
already in Refs. [4, 27]. All terms given above agree nicely with our previous group 
theoretical analysis. 



2.7. Numerical estimates 

We now estimate the size of the Rashba spin-orbit coupling terms in silicene. The 
bond parameters Vi_3, the energy d, the angle 9 = 101.7°, the lattice constant a, as 
well as the spin orbit interaction A are taken from Ref. [12]. For the Rashba-terms, 
we then calculate the Stark-element zq = {(f)n,ofi\z\(f)n,i,o) as transition matrix elements 
between atomic wave functions, where in silicene we have n = 3, since the outer shells 
are provided by the 3s and 3p orbitals. The corresponding wave functions were chosen 
as the wave functions of the hydrogen atom with the same quantum numbers. We then 
find Zq = 3\/6 x aB/Zes in silicene, where is the Bohr radius and an outer electron 
experiences an effective atomic charge of Z^s ~ 4.29 [28] due to screening. Hence, we 
estimate Zq = 0.906A. For typical values of E^ = 50eV/300nm, we calculate 

Ae ^ 8.5 meV, (23) 
Ae,2 ~ 12.8 ficV, 
Xr,i ^ 22.7 fieV. 

Thus, our new term proportional to Ae,2 is small compared to the term proportional to 
Ae, but similar in magnitude as the Kane-Mele-Rashba term proportional to Ar,i. 



3. Topological analysis 

Triggered by the theoretical prediction [4, 5, 6] and experimental discovery [8] of the 
quantum spin Hall (QSH) effect, the study of topological effects in the physics of Bloch 
bands has been a major focus of condensed matter physics in recent years [1, 2, 3]. The 
QSH state is a bulk insulating state featuring metallic edge states that are protected 
by time reversal symmetry (TRS). Due to Kramers theorem, these edge states appear 
in so called helical pairs. The bulk energy gap in the original proposal for the QSH 
effect in graphene by Kane and Mele [4] is due to an intrinsic SOI which preserves 
a residual U{1) spin symmetry. In the presence of this global spin quantization axis. 
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the helical edge states are characterized by a perfect locking of spin and momentum: 
States with opposite spin move with opposite chirality along the edge. In the presence 
of Rashba SOI no spin symmetry is present resulting in the absence of a global spin 
quantization axis. Rather, the spin quantization axis of the edge states can precess 
spatially. However, as TRS is still present the Kramers pair of helical edge states is 
still protected from hybridizing, i.e., from a gap opening on the edge. This statement 
is true as long as the bulk gap is maintained implying that the system is adiabatically 
connected to the U{1) preserving case and hence is still in the QSH state. However, 
the role of Rashba SOI for the QSH effect in graphene is only potentially detrimental: 
Upon increasing the strength of Rashba SOI, the system will go through a quantum 
phase transition destroying the QSH phase. Conversely, given the lattice symmetry 
of graphene, switching on Rashba spin orbit coupling cannot generate the QSH phase 
starting from a semi-metallic or trivial insulator phase. 

The slightly reduced symmetry of silicene stemming from its buckled structure 
entails several new couplings, at least two of which are of key relevance regarding QSH 
physics. First, let us consider the staggered potential distinguishing sub lattice A and 
B that has been introduced formally in Ref. [4] to open a trivial gap. While this term 
is symmetry forbidden in graphene, it has been shown [14] to be induced by a simple 
out of plane electric field in silicene as we confirmed in our symmetry analysis above. 
This provides a knob to experimentally switch off the QSHE. Second, and even more 
interestingly, it has very recently been conjectured [15] that the QSHE can even be 
generated by virtue of Rashba SOI terms that are symmetry forbidden in graphene 
but allowed in silicene. The authors of Ref. [15] probe the sub-gap conductance as 
a fingerprint of the quantum spin Hall state. While a non-vanishing conductance in 
a phase with a bulk gap is a promising signature of the QSH state it is not in one 
to one correspondence with the Z2-invariant defining the quantum spin Hall effect [5]. 
For example, what is called the QSHE2-phase in Ref. [15] is a trivial insulator which 

2 

only shows the characteristic conductance of 2^ expected for the QSH phase since one 
additional pair of edge states does not contribute due to a mini-gap at the ribbon sizes 
considered in Ref. [15]. In the thermodynamic limit, the conductance in this parameter 

2 

regime would be 4y signaling a Z2-trivial phase. 

The main purpose of this section is to calculate the correct phase diagram of silicene 
in the presence of the Rashba SOI terms by rigorous calculation of the Z2 invariant in 
the absence of any symmetries besides TRS. We note that our method can be applied 
to study the entire parameter space of the silicene band structure without any further 
complications. However, in this work, we would like to focus on a particularly interesting 
parameter regime where the QSH phase is generated by Rashba SOI. Previous work on 
the topology of the QSH state in silicene in the presence of Rashba SOI [14] has been 
focused on effective models vahd close to the K-points. However, in the presence of 
Rashba SOI, the bulk gap can close away from the K-points, a phase transition which 
would be missed by such power series expansions. Our analysis does not suffer from 
these limitations since we directly calculate the Z2-invariant characterizing the QSH 
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Hall phase from its very definition [5, 17]. We find that there is indeed a QSH phase 
in the absence of the original Kane-Mele term A^o- This QSH effect can be switched on 
by only tuning silicene specific Rashba terms and the staggered potential-all terms that 
are known to have only detrimental effect as to QSH physics in graphene. Our analysis 
hence settles the discussion whether Rashba SOI can in principle generate a QSH state 
in the affirmative. 

3.1. Manifestly gauge invariant calculation of the topological 7j2-invariant 

Let us briefly give an idea of how the topology of silicene may depend on an external 
electric field, orienting ourselves along the lines of Ref. [14]. We consider again the 
effective Hamiltonian of Eq. (11), derived by analytical expansion around the high- 
symmetry i^-points. Without electric fields, the energy spectrum at K exhibits a gap 
of size |2a2|. Interestingly, in the presence of SOI, the bulk gap can be closed by an 
increasing electric field perpendicular to the plane. The gap closes approximately at 
the very i^-points, only if the corresponding parameter A^,! is small compared to the 
transfer energy t, so A^.i << t. The low-energy Hamiltonian allows to give an analytical 
expression for the gap-closing critical field [14]. At this point we may expect a 
topological phase transition from a QSH to a trivial insulating phase. This is indicated 
by the fact, that a gap is reestablished at the i^-points, if the electric field is further 
increased, exceeding E';. However, a rigorous exploration of prospective topological 
quantum phase transitions in the silicene parameter space requires the calculation of a 
topological bulk invariant. 

Therefore, we now turn to the general calculation of the topological Z2-invariant 
defining the QSH phase [5]. While the Z2-invariant is of course a gauge invariant 
quantity by definition, the original literature [5, 16] did not provide a constructive 
recipe for its numerical calculation. This is because a calculation following the original 
definition requires a macroscopic gauge, though giving a gauge invariant result. By 
macroscopic gauge, we mean that the phase relation between Bloch functions at remote 
points in the Brillouin zone has to be fixed. However, when the band structure is 
calculated numerically, such phase relations are typically not accessible thus preventing 
the direct calculation of the invariant. This problem has only been resolved rather 
recently [17, 29] by a more constructive recipe for the direct calculation of the Z2- 
invariant. Here, we follow the method introduced by Prodan [17] which makes use 
of the elegant and manifestly gauge invariant formulation of the adiabatic connection 
[3, 30, 31, 32] originally introduced in a seminal work by Kato back in 1950 [30]. Recently, 
the manifestly gauge invariant formulation of topological invariants has been generalized 
to other topological band structures [3]. The crucial step of this approach consists 
in going from the Bloch functions of the occupied bands to the projection operator 
P{k) onto the occupied states. As already pointed out by Kato in Ref. [30], the 
advantage of this construction is that the projection operator is obviously basis (gauge) 
independent. 
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The adiabatic connection is defined as 

A,{k) = -[id,P{k)),P{k)], (24) 

where 9^ = is the derivative in momentum space. Note that this connection is 
manifestly independent of the basis choice within the occupied bands in contrast to the 
more famihar Berry connection. Therefore, the adiabatic connection can be calculated 
numerically in a straight forward way which is in general not possible for the Berry 
connection. In Ref. [17], the author constructs the Z2 invariant of the QSH state in 
a similar way to Ref. [16] but using the adiabatic connection instead of the Berry 
connection. We refer the reader to Ref. [17] for this both accessible and fairly self 
contained explicit construction and review only the resulting expression for the Z2- 
invariant H for a rectangular lattice with lattice constants a^, ay, for completeness, 

^ ^ ^ det(?7(o,o),(o,bj) det(?7(j,^,o),(fc^,b^))Pf(g(o,o))Pf(g(j,^,o)) ^^^^ 
Pf(^(o,fe,))Pf(^(b.,fe,))Vdet(f/(o,b^),(o,-b,))Vdet(f/(b^,b^),(b^,_f,^)) 

Here, = f^? ~ ^ boundaries of the Brillouin zone, ^(fc^,/cj,) is the matrix 

of the TRS operation which is anti-symmetric at the time-reversal invariant momenta 
(TRIM), and Pf denotes the Pfaffian. We note that S = — 1 defines the non-trivial QSH 
phase whereas S = 1 for a trivial insulator, a will be explained in detail below and 

f^(;tw,fc(/)) = re-4w 

describes the adiabatic evolution along the straight line from /c*^*^ to k^^^ with the path 
ordering operator T. The gauge invariance of Eq. (25) might not be obvious at first 
glance due to the basis dependence of the Pfaffians. However, the combinations 



f det(?7(fc„o),(fc,,b,))Pf(g(fc„o)) f—— : 

= ^77^ ^ = ±y <iet[U(^k^^by),ik^,-by)) 1^6) 

with kx = 0, bx appearing in Eq. (25) are indeed gauge invariant up to the choice of the 
branch of the multivalued square-root as is obvious from the right hand side of Eq. (26). 
That is the point where a comes into play. When calculating 

^/det{U^o,by),(o-by))^/det{U(^bx,by),{bx-by) 
it has to be assured that the branch choices of the square-root at fc^, = and k^ = b^ are 
continuously connected to each other. This can be done by continuously interpolating 
the phase factor det{U(^kx,by),(kx-by)) between k^ = and k^ = b^ [17]. If this phase has 
a winding such that the standard branch cut between Riemann sheets (line (—00, 0] in 
the complex plane) is crossed n times during this interpolation, the naive result for S is 
corrected by the factor a = (—1)"". 

While this procedure might be numerically challenging close to critical points or 
for large super-cells in disordered systems, it is basically a straightforward recipe. Note 
that the arbitrary basis choice for the representation matrix Oi^k^My) does not require any 
knowledge about the relative phase between the basis vectors at different points in k- 
space. Furthermore, the mentioned phase interpolation does not require any numerically 
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inaccessible information either, since the phase factor det{U(^k^,by).{k^,-by)) at every is 
a gauge invariant quantity. 

3.2. From the honeycomb lattice of silicene to a rectangular super-cell 

The procedure for the calculation of the Z2-invariant just described is general but 
requires a rectangular lattice whereas silicene crystalizes in a honeycomb lattice. On the 
other hand, the result for the topological invariant cannot depend on the choice of the 
unit cell. Therefore, we will now introduce a rectangular super cell which immediately 
allows us to use the above method. To this end, the common basis vectors given in 
Eq. (1) are combined to a new set of basis vectors 

d[ = di + 0,2 = (0, 3a, 0) , a'g = ai — 02 = (^VSa, 0, , (27) 

spanning a rectangular lattice with four atoms per site A, B, A' and B', as shown in 
Fig. 3. The Brioullin zone is again rectangular with basis vectors 

S;.(o.|,o). 9.-{^/>,o). (28) 





Figure 3. Unit cells of the two-dimensional honeycomb lattice. On the left hand 
side, the minimal unit cell is drawn, containing two basis atoms A and B. On the right 
hand side, the unit cell has a rectangular form of doubled size. Here, four basis atoms 
A,B, A' , and B' are included. 

The Bloch Hamiltonian of size (8 x 8) associated with our tight-binding model 
consists of the following contributions 

H = Ht + H,o + H, + H,,2 + Hx^, + H^^, . (29) 

The form of the fc-dependent matrices is listed in Eqs. (A.l) to (A. 6) in Appendix A. 
In the extended unit cell, the i^-points (former corner points of the Brillouin zone) are 
mapped onto points (±^^^^,0) inside the rectangle forming the new Brillouin zone. 

For the Hamiltonian we chose the basis {IV^a), H^b), li^A'), iV's')} ^ {t?!}; where the 
wavefunctions are of the form [33] 
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with X G {A, B, A', B'}. cf) are atomic wavef unctions, in particular {r\(l)f) = (f){f— Rf)- 
Rj denote lattice vectors in real space, connecting the reference points of unit cells, while 
Rf are the positions of atoms labeled X, in the unit cell j, relative to the reference 
point. The sum runs over all cells of the crystal. 

3.3. Results for the phase diagram of silicene 

For notational simplicity, the electric field is absorbed in the parameters Aj, that is 
XiEz — )■ Aj. All parameters are measured in units of the hopping t. We keep in mind 
that Ag, Ae,2, and A^,! depend on the external electric field. 

We first benchmark the general method in a parameter regime where we have a 
good intuition for the expected QSH physics from previous literature [4, 14]. Let us 
start with just one non-vanishing parameter t, all other parameters being zero. This 
produces the well known gapless Dirac cones at the i^-points. Upon switching on A^o 
[4], the system exhibits a bulk gap of size \2Xso\ and we reproduce S = —1 [5]. We now 
add a term proportional to Ag depending on the external field. Upon increasing Ag the 
bulk gap closes at the /^-points for a critical value as predicted before. For even stronger 
fields a trivial gap opens, i.e., S = +1. Similarly, the gap closes for terms proportional 
to Ae,2- Next, we add the Rashba-terms proportional to Ar,i and Ar,2- The first term Ar,i 
tends to establish three additional minima of the energy gap, that are shifted away from 
the i^'-point. On the other hand, Ar._2, primarily providing gapless states at the corners 
of the BZ, causes a modification of the band slope. From an analysis of the overall band 
structure, one finds, that the interplay of both terms Ar.,i and Ar,2 may possibly close 
the band gap away from the if-points. After these general considerations, we would 
like to focus on the interesting parameter regime identified in Ref. [15], i.e., we consider 
Ae, Xr,i, Ar,2 as free parameters that are measured in units of t. All other parameters 
are set to zero. Unfortunately, the analysis in Ref. [15] was not suitable to predict 
the correct phase diagram for the Z2 invariant. However, our more rigorous analysis 
confirms the general conjecture that in this regime, a combination of A^,! and Ar,2 is able 
to drive a topological quantum phase transition away from the if-points which induces 
a non-trivial QSH phase. If we only switch on Ag, a trivial gap of size |2Ag| opens. 
By tuning Ar,i and Ar,2, the bulk gap can be closed away from the i^'-points, and a gap 
characterized by S = —1, i.e., a QSH state emerges. We present the phase diagram in 
the identical parameter regime as presented in Ref. [15] in Fig. 4. Our direct calculation 
of the Z2-invariant disagrees with Ref. [15] both qualitatively (absence of the QSHE2) 
and quantitatively (significantly different phase boundary for the QSH phase). However, 
we would again like to point out that we can definitely confirm the phenomenology of a 
QSH phase in the absence of the Kane-Mele term Aso- Instead, the extended QSH region 
shown in Fig. 4 is only driven by the two silicene specific Rashba SOI terms A^,! and 
Ar,2 which is conceptually very interesting. 

Lastly, we found the phase diagram in Fig. 4 to be robust against small perturbations 
proportional to the term A^q. For larger Kane-Mele parameters A^o > Ag/ (3v^), another 
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non-trivial regime (S = —1) enters the phase diagram. This regime occurs independently 
of Ar,2 and can be abandoned again by increasing Xr^i- Essentially, the previously 
described quantum phase transitions are not affected. 




0.4. 0.6 



1.0 



Figure 4. Topological phase diagram of silicene with free parameters A^.i, Xr.2- The 
staggered potential Ag = 0.1 is fixed as well as Xgo = 0. All couplings are measured in 
units of the hopping energy t. The white points close to the phase boundary denote 
critical regions where the bulk gap is so small that our numerical calculation did not 
converge properly. 



4. Conclusion 

A symmetry analysis of the lattice of silicene close to the K points of the Brillouin zone 
was performed. With the help of the invariant expansion model, the vr-band Hamiltonian 
was constructed by symmetry considerations only, including spin-orbit coupling and 
external electric fields perpendicular to the atomic plane. Thereby, we discovered a new 
term that is allowed by symmetry and related to an interplay of spin-orbit coupling 
and an external electric field in a buckled honeycomb structure. We supplemented the 
symmetry analysis by a tight-binding model which allowed us to estimate the relevant 
parameters of the model. 

Subsequently, we carefully analyzed the topological properties of the band structure 
and proved that a topological phase transition can be generated by an external 
electric field. This analysis enabled us to plot a topological phase diagram of silicene 
employing Prodan's manifestly gauge-invariant method for the direct calculation of the 
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Z2 topological invariant defining the quantum spin Hall phase. Since this method 
requires a rectangular lattice, we considered a rectangular silicene super-cell that 
contains four atoms. Interestingly, the tunable phase transition can happen in a 
destructive as well as a constructive way, i.e., a quantum spin Hall phase can not only 
be destroyed but also be generated by means of external Rashba spin orbit coupling. 
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Appendix A. Tight-binding model with rectangular unit cell 



The terms contributing to the tight-binding Hamiltonian in Eq. (29), in the basis 
{\iPa), \^Pb),\'^A'),\'^b')} X {t,i}, take the explicit form 
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The matrix elements are given in tables Al and A2. 
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Table Al. Matrix elements in Eqs. (A.l) to (A. 6) of the diagonal matrix blocks. 
We use the abbreviations x = ^^k^ and y = ^ky. 

2t cos(3;) (cos{3:) + z sin{a:) ) ui ^ lA^ 2 cos(x) sin(3::) 

-2Aao ain(2a;) ei X^E^ 

(cos(a;) + 2 sin(a;)) (cos(x) + \/3sin(3^)) 



\/l + cot2(9) 



yi+cot2(£)) 

(cos(a;) + I sin(a;))(cos(ic) — \/3sin(3:;)) 



Table A2. Matrix elements in Eqs. (A.l) to (A. 6) of the non-diagonal matrix blocks. 

^^-^kx and y = ^ky. 



We use the abbreviations x = ^^^^k^ and y = ^ky 



t(cos(6y) + isin.{Qy)) t.-, t 

" -^e 2 

4Aso sin(x) cos(3y) [cos (x + 3y) + i sin (x + 3y)] £2 ■^i^l — ^5; — 

/ 2 ^ so e 

"^2^1 A — A — ^2 4Aso sin(x) cos(3y) [— cos (a; — 3y) + i sin (x — 3y)] 



^ ^ = (sin{6y) — icos(6t/)) 



^l + cot2(0) y/l + cot2(e)} 

(i (V3 + 3i) sin(x - 3t/) + (3 + iVs) sin(a:: + 31/)) X ((^3 - 3i) sin(cr - 3y) + (V3 + 3i) sin(ic + 3y)) X 

(cos(a:: + 3y) + z sin(x + 3y)) (sin(x + 3y) — i cos(x + 3y)) 

((3 - zn/s) sin(x - 3y) + (-3 - zVs) sin(x + 31/)) X ((%/3 - 3z) sin(x - 3y) + (^3 + 3i.) sin(a; + 3y)) X 

(cos(a: — 3y) — z sin(a: — 3y)) (sin(x — 3y) + z cos(x — 3y)) 



Note, that coefficients are related to parameters of the group theoretical approach by 
Oi = 3a/3A<jo, 02 = A(= 



^3 



as 



3\/3Ae2) (^4 ~ Af 1 



2v/l + cot2(0)' 
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